/*
 * Dipoles.cpp
 *
 *  Created on: Sep 4, 2013
 *      Author: marchi
 */

#include "Dipoles.h"
vector<string> split(const string & s){
	stringstream iss(s);
	vector<string> tokens;
	copy(std::istream_iterator<string>(iss),
			std::istream_iterator<string>(),
			std::back_inserter<vector<string> >(tokens));
	for(unsigned int o=0;o<tokens.size();o++)
		tokens[o].erase(remove_if(tokens[o].begin(),tokens[o].end(),::isspace),tokens[o].end());
	return tokens;
}


void Dipoles::initLists(TopolPDB & y,vector<string> & L){
	map<string,map<string,double> >Charges;
	this->MAtoms::initLists(y,L);
	vector<string> data;
	for(string str;getline(*fin,str);){
		data.push_back(str);
	}
	size_t o=0;
	do {
		vector<string> title=split(data[o]);
		int found=title[0].substr(0,1).find("$");
		if(found != string::npos){
			o++;
			map<string,double> Res;
			vector<string> tokens=split(data[o]);
			do{
				stringstream ss(tokens[1]);
				double tmp;
				ss>>tmp;
				Res[tokens[0]]=tmp;
				o++;
				tokens.clear();
				tokens=split(data[o]);
			} while(tokens[0].substr(0,1).find("$") == string::npos  );
			string sub1=title[0].substr(1,title[0].size()-1);
			Charges[sub1]=Res;
		} else{o++;}
	} while (o<data.size() && data[o].find("$end") == string::npos);
	chg=vector<double>(nr,0.0);
	for(size_t o=0;o<PDB.size();o++){
		string atn=PDB[o].atn;
		string resn=PDB[o].resn;
		if(Charges.count(resn) > 0){
			if(Charges[resn].count(atn) > 0){
				chg[o]=Charges[resn][atn];
			}
		}
	}
}
void Dipoles::gDipole(){
	Dvect dip=0.0;
	for(size_t o=0;o<ResIndx0->size();o++){
		size_t i=(*ResIndx0)[o];
		Dvect x0=x[i];
		x0*=chg[i];
		dip+=x0;
	}
	dipole.push_back(dip);
}
void Dipoles::dPrint(ostream & fout){
	double fact1=4.0*M_PI*(Parameters::efact/1000.0/Parameters::kT300)/Parameters::unitc;
	cout << fact1 <<endl;
	vector<double> Mycorr=vector<double>(dipole.size(),0.0);
	vector<int> iMycorr=vector<int>(dipole.size(),0);
	for(size_t m=0;m<dipole.size();m++){
		for(size_t o=0;o<dipole.size()-m;o++){
			Mycorr[m]+=dipole[o]*dipole[o+m];
			iMycorr[m]++;
		}
	}
	/*
	Gyration Rg=Rg_stat;
	Rg=Rg/Rg_count;
	cout << Rg.Volume() << " " << Rg.Axis()[0] << " " << Rg_count << endl;
	double Vol=Rg.Volume()/1000.0;
	for(size_t m=0;m<dipole.size();m++){
		double fact=static_cast<double> (iMycorr[m]);
		Mycorr[m]/=fact;
	}
	double tmp0;
	tmp0=Mycorr[0]*fact1/Vol;
	double eps=1.0+tmp0;
	fout << " Eps = " << eps << "  " <<endl;
	cout << " Eps = " << eps << "  " <<endl;
	for(size_t m=0;m<dipole.size();m++){
		fout << m << " " << Mycorr[m] <<endl;
	}
	*/
}
Dipoles::~Dipoles() {
	// TODO Auto-generated destructor stub
}

